######################################################################
###################  	 REPLICATION CODE	  ########################


# Lucas Leemann, 3.11.2015
# l.leemann@ucl.ac.uk



rm(list=ls())
require(MCMCpack)
require(MSBVAR)
require(xtable)
require(foreign)
require(MASS)
library(memisc)
library(caTools)
library(pingr)
library(texreg)

my_path <- "write_your_path_here"
setwd(my_path)
mydata    <-  read.dta("Data_SPSR_2015_Leemann.dta")


mydata[c(90:93),2] <- c(5,5,1,6) 		# years 2008, 2009, 2010, 2011
								# added from https://www.admin.ch/ch/d/pore/vr/vor_2_2_6_3_13_2011.html
mydata[c(90:93),23] <- c(4,5,3,2) 		# years 2008, 2009, 2010, 2011
								# added from http://www.admin.ch/ch/d/pore/vr/vor_2_2_6_3_13_2001.html

######################################################################
###################### 		CODEBOOK		##########################

# year 
#	> Year

# number_of_DD 
#	> Number of Initiatives SUBMITTED in this year (swissvotes.ch and admin.ch)

# election 
#	> Binary indicator whether it was an election year (=1) or not (=0)

# change_PartyCanton 
#	> Absolute change in party support in %-points within a canton. Value is mean over all cantons and refers to change
#	> from last elections to next election (refers to period over which the change occurs). Value in year *t* is based on last election *t-x* (where x is maximally 3) and
#	> *t+y* (where y is maximally 3). Maximal value is 100.

# change_partyNation
#	> Same as change_PartyCanton, but calculated based on national vote returns


# _merge merL eme mergI
#	> merging variables

# Year2 Year3
# 	> Year squared and year cubed.

# post1970 post1990
#	> Indicator variables indicating whether a year is after 1970 (1990)

# Ychange
#	> Interaction between Year and change_PartyCanton

# right_to_vote 
#	> number eligible voters, source: Bundesamt fÃ¼r Statistik, Statistik der eidg. Volksabstimmungen. Â© BFS - Statistisches 
#	> Lexikon der Schweiz http://www.bfs.admin.ch/bfs/portal/de/index/themen/17/03/blank/key/stimmbeteiligung.html 

# requirement 
#	> signature requirement, http://www.hls-dhs-dss.ch/textes/d/D10386.php

# share 
#	> Signature requirement divided by vote entitled citizens (own calculations), in % points  

# legislature 
#	> Legislative period of the National Council

# seatshare_change 
#	> Change in seat shares by parties in the National Council in % (max value 1)

# time
#	> year - 1915 = time

# pull_back
#	> number of initiatives which were submitted in that year but eventually pulled back.




######################################################################
#			data manipulation
######################################################################

# getting rid of 1919 - will not be covered in the analysis
mydata <- mydata[-1,]


#normalizing some variables to avoid convergence problems later
change_pN <- (mydata$change_partyNation-mean(mydata$change_partyNation))/sd(mydata$change_partyNation)
change_pC <- (mydata$change_PartyCanton-mean(mydata$change_PartyCanton))/sd(mydata$change_PartyCanton)
sharE  <- (mydata$share-mean(mydata$share))/sd(mydata$share)
SS_change <- (mydata$seatshare_change-mean(mydata$seatshare_change))/sd(mydata$seatshare_change)
loss_pC <- (mydata$loss_PartyCanton-mean(mydata$loss_PartyCanton))/sd(mydata$loss_PartyCanton)
part_st <- (mydata$participation-mean(mydata$participation))/sd(mydata$participation)

# create a moving average
ma.DD <- rep(NA, length(mydata$number_of_DD))
ma.DD[1] <- mean(mydata$number_of_DD[c(1:2)])

for (i in 2:92){
		c <- i-1
		d <- i+1
		ma.DD[i] <- mean(mydata$number_of_DD[c(c:d)])
}
ma.DD[92] <- mean(mydata$number_of_DD[c(91:92)])


# number of atcual initiatives submitted (subtracting those that were pulled back -- used in robustness section)
number_of_DD1 <- mydata$number_of_DD + mydata$pull_back



############################################################################################################################################
#			Empty Count models 
############################################################################################################################################

mcmc <-  200000
burnin <-  100000
thin   <-  100
verbose <-  10000
out.empty0 <-  MCMCpoisson(number_of_DD1 ~ 1, data=mydata, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin, b0=0, B0=0.01, marginal.likelihood = "Laplace")  

burnin <-  200000                           
out.empty1 <-  MCMCpoissonChange(number_of_DD1 ~ 1, data=mydata, m=1, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin,c0=4.6, d0=1, marginal.likelihood = "Chib95")  

mcmc <- 300000
out.empty2 <-  MCMCpoissonChange(number_of_DD1 ~ 1, data=mydata, m=2, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin,c0=3.1, d0=1, marginal.likelihood = "Chib95")  
                           
mcmc <- 600000
out.empty2 <-  MCMCpoissonChange(number_of_DD1 ~ 1, data=mydata, m=2, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin,c0=3.1, d0=1, marginal.likelihood = "Chib95")                        

burnin <-  1500000
mcmc <- 1000000
out.empty3 <-  MCMCpoissonChange(number_of_DD1 ~ 1, data=mydata, m=3, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin,c0=2.3, d0=1, marginal.likelihood = "Chib95")  
burnin <-  100000
mcmc <-  200000
out.empty4 <-  MCMCpoissonChange(number_of_DD1 ~ 1, data=mydata, m=4, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin,c0=1.4, d0=1, marginal.likelihood = "Chib95")  

# checking convergence
geweke.diag(out.empty0)
geweke.diag(out.empty1)
geweke.diag(out.empty2)
geweke.diag(out.empty3)
geweke.diag(out.empty4)


heidel.diag(out.empty0)
heidel.diag(out.empty1)
heidel.diag(out.empty2)
heidel.diag(out.empty3)
heidel.diag(out.empty4)


# Regimes
par(mfrow=c(2,1))
plotState(out.empty2)
plotState(out.empty3)

# Bayes Factor to determine support by data
print(BayesFactor(out.empty0, out.empty1, out.empty2, out.empty3, out.empty4))
log.bayes <- BayesFactor(out.empty0, out.empty1, out.empty2, out.empty3, out.empty4)$BF.log.mat
colnames(log.bayes) <- c("No Break", "One Break", "Two Breaks", "Three Breaks", "Four Breaks")
rownames(log.bayes) <- c("No Break", "One Break", "Two Breaks", "Three Breaks", "Four Breaks")
xtable(log.bayes)


prob.state <- attr(out.empty2, "prob.state")    
change <- prob.state[,2]
change1 <- abs(diff(prob.state[,1]))# + abs(diff(prob.state[,2]))
change2 <- abs(diff(prob.state[,2]))# + abs(diff(prob.state[,3]))

Schange <- c(change1[1:40],change2[41:91])


############################################################################################################################################
#			F	I	G	U	R	E 		1
############################################################################################################################################

# Best model
mcmc <- 700000
out.empty2 <-  MCMCpoissonChange(number_of_DD1 ~ 1, data=mydata, m=2, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin,c0=3.1, d0=1, marginal.likelihood = "Chib95")                        



###############################################################################################################
# FIGURE 1 (color)
par(family="CMU Serif",mai=c(0.9,1,.1,.3),oma=c(.4,0,.1,0))
layout(c(1,2),heights=c(3,1))
plot(mydata$year, mydata$share, ylim=c(0,0.10), xlim=c(1919,2012), type="l", lwd=4, lty=2, yaxt="none", xlab="", bty="n", col="white", xaxt="none", ylab="")
rect(1919,0.00,2012,0.1015, col=rgb(247,247,247,150,max=255), border=NA)

points(mydata$year, number_of_DD1/100,col=rgb(30,144,255,150,max=255),cex=0.7,pch=19)
	for (i in seq(100,600,length.out=150)){
		points(lowess(mydata$year, jitter(number_of_DD1, factor=2)/100,f=i/1000), type="l",col=rgb(30,144,255,20,max=255), lwd=.99)
	}

axis(2, at = c(0,0.02,0.04,0.06,0.08, 0.1), labels=c("0","2","4","6","8","10"), col.lab=rgb(30,144,255,250,max=255), col.axis=rgb(30,144,255,250,max=255) , col=rgb(30,144,255,250,max=255), las=1, cex.axis=1.3)

points(mydata$year[80], 0.1,col=rgb(30,144,255,150,max=255),cex=1.0,pch=17)
axis(1,at=c(1919,1930,1940,1950,1960,1970,1980,1990,2000,2010), labels=c("1921","1930","1940","1950","1960","1970","1980","1990","2000","2010"), line=.1,xpd=NA,cex.axis=1.2)#pos=c(-.1,0))
mtext(side=4,"Number of Submitted Initiatives", line=-.8, col=rgb(30,144,255,200,max=255), cex=1.6)


#par(family="CMU Serif", mai=c(0.1,1,0.8,.3))
barplot(Schange/1.8, yaxt="n", width=rep(.5,91), border=NA, col=rgb(255,64,64,200,max=255), bg=rgb(247,247,247,150,max=255))
rect(0,0.00,55,0.32, col=rgb(247,247,247,150,max=255), border=NA, xpd=NA)
barplot(Schange/1.8, yaxt="n", width=rep(.5,91), border=NA, col=rgb(255,64,64,200,max=255), bg=rgb(247,247,247,150,max=255), add=TRUE)

axis(2,at=c(0,0.06,0.11,0.17,0.22,0.27,.32), labels=c("0%","","20%","","40%","","60%"), xpd=NA, cex.axis=1.1, col.lab=rgb(255,64,64,200,max=255),col.axis=rgb(255,64,64,200,max=255),col=rgb(255,64,64,200,max=255),las=1)

mtext(side=4,"     Probability of \n      Regime Change", line=-.18, col=rgb(255,64,64,200,max=255), cex=1.5)


###############################################################################################################
# FIGURE 1 (black and white)
par(family="CMU Serif",mai=c(0.9,1,.1,.3),oma=c(.4,0,.1,0))
layout(c(1,2),heights=c(3,1))
plot(mydata$year, mydata$share, ylim=c(0,0.10), xlim=c(1919,2012), type="l", lwd=4, lty=2, yaxt="none", xlab="", bty="n", col="white", xaxt="none", ylab="")
rect(1919,0.00,2012,0.1015, col=rgb(247,247,247,150,max=255), border=NA)

points(mydata$year, number_of_DD1/100,col="grey",cex=0.7,pch=19)
	for (i in seq(100,600,length.out=150)){
		points(lowess(mydata$year, jitter(number_of_DD1, factor=2)/100,f=i/1000), type="l",col=rgb(150,150,150,40,max=255), lwd=.99)
	}

axis(2, at = c(0,0.02,0.04,0.06,0.08, 0.1), labels=c("0","2","4","6","8","10"),  las=1, cex.axis=1.3)

points(mydata$year[80], 0.1,col=rgb(190,190,190,150,max=255),cex=1.0,pch=17)
axis(1,at=c(1920,1930,1940,1950,1960,1970,1980,1990,2000,2010), labels=c("1920","1930","1940","1950","1960","1970","1980","1990","2000","2010"), line=.1,xpd=NA,cex.axis=1.2)#pos=c(-.1,0))
mtext(side=4,"Number of Submitted Initiatives", line=-.8, cex=1.6)


#par(family="CMU Serif", mai=c(0.1,1,0.8,.3))
barplot(Schange/1.8, yaxt="n", width=rep(.5,91), border=NA, col=rgb(150,150,150,230,max=255), bg=rgb(247,247,247,150,max=255))
rect(0,0.00,55,0.32, col=rgb(247,247,247,150,max=255), border=NA, xpd=NA)
barplot(Schange/1.8, yaxt="n", width=rep(.5,91), border=NA, col=rgb(150,150,150,230,max=255), bg=rgb(247,247,247,150,max=255), add=TRUE)

axis(2,at=c(0,0.06,0.11,0.17,0.22,0.27,.32), labels=c("0%","","20%","","40%","","60%"), xpd=NA, cex.axis=1.1, las=1)

mtext(side=4,"     Probability of \n      Regime Change", line=-.18, cex=1.5)






############################################################################################################################################
#			Count models with no change
############################################################################################################################################

options(signif.symbols=c("***"=.01,"**"=.05,"*"=.1))

# generating pre-election variable
pre.election <- rep(NA,length(mydata$election))
pre.election[1:91] <- mydata$election[2:92]
pre.election[92] <- 0

# generating post-election variable
#post.election <- rep(NA,length(mydata$election))
#post.election[2:92] <- mydata$election[1:91]
#post.election[1] <- 1

# Table 1 paper
formula <- number_of_DD1 ~ change_pC + sharE + SS_change + election + pre.election 

nbr1 <-  glm.nb(number_of_DD1 ~  change_pC , data = mydata)
nbr2 <-  glm.nb(number_of_DD1 ~  sharE  , data = mydata)
nbr3 <-  glm.nb(number_of_DD1 ~  change_pC + sharE  , data = mydata)
nbr4     <-  glm.nb(number_of_DD1 ~   change_pC + sharE + election , data = mydata)
nbr5     <-  glm.nb(number_of_DD1 ~   change_pC + sharE + election + pre.election , data = mydata)

table1 <-   mtable("Model 1"=nbr1,"Model 2" = nbr2, "Model 3" = nbr3, "Model 4" = nbr4, "Model 5" = nbr5, digits=2, summary.stats=c("McFadden R-sq.","N","Likelihood-ratio","p","BIC")) 
table1 <- relabel(table1, "(Intercept)"="Constant", "change_pC"="Electoral Competition", "sharE"="Signature Requirements", "election"="Election Year", "pre.election"="Pre-Election Year")
toLatex(table1)

# Robsutness model 1.1

nbr1.r <-  glm.nb(number_of_DD1 ~  change_pN , data = mydata)
nbr2.r <-  glm.nb(number_of_DD1 ~  sharE  , data = mydata)
nbr3.r <-  glm.nb(number_of_DD1 ~  change_pN + sharE  , data = mydata)
nbr4.r     <-  glm.nb(number_of_DD1 ~   change_pN + sharE + election , data = mydata)
nbr5.r     <-  glm.nb(number_of_DD1 ~   change_pN + sharE + election + pre.election , data = mydata)

table1 <-   mtable("Model 1"=nbr1.r,"Model 2" = nbr2.r, "Model 3" = nbr3.r, "Model 4" = nbr4.r, "Model 5" = nbr5.r, digits=2, summary.stats=c("McFadden R-sq.","N","Likelihood-ratio","p","BIC")) 
table1 <- relabel(table1, "(Intercept)"="Constant", "change_pN"="Electoral Competition (national)", "sharE"="Signature Requirements", "election"="Election Year", "pre.election"="Pre-Election Year")
toLatex(table1)
screenreg(list(nbr1.r,nbr2.r,nbr3.r,nbr4.r,nbr5.r),include.pvalues=TRUE ,stars = c(0.01,0.05,0.1))




############################################################################################################################################
#			F	I	G	U	R	E 		2
############################################################################################################################################


############################################################################################################################################
# color
dataDim    <-  read.dta("Data_SPSR_2015_Leemann_dim.dta")
dataDim <- dataDim[-dim(dataDim)[1],]
dataDim1 <- rbind(dataDim,c(2011,0,2,0,0,1,1,0,0))
## Plot
par(family="CMU Serif", mai=c(1.0,1,0.5,1.0))#,bg="grey95")
plot(dataDim1$year,dataDim1$Dimension1, ylim=c(0,4.9), xlim=c(1919,2010), type="l", lwd=2, lty=2, yaxt="none", xlab="Year", bty="n", col=rgb(255,64,64,1,max=255), xaxt="none", ylab="Number of Annual Initiatives")
rect(1919,-0.1,2012,6.4, col=rgb(247,247,247,150,max=255), border=NA)
points(jitter(dataDim1$year),dataDim1$Dimension2-0.04, col=rgb(255,64,64,50,max=255), pch=19)
points(jitter(dataDim1$year),dataDim1$Dimension1+0.04, col=rgb(30,144,255,50,max=255), pch=19)
#
points(lowess(dataDim1$year,jitter(dataDim1$Dimension1), f=0.2), col=rgb(30,144,255,200,max=255), type="l", lwd=4, lty=4)
for (t in 120:350){
	points(lowess(dataDim1$year,jitter(dataDim1$Dimension1), f=t/1000), col=rgb(30,144,255,20,max=255), type="l", lwd=.4, lty=1)	
}
points(lowess(dataDim1$year,jitter(dataDim1$Dimension2), f=0.2), col=rgb(255,64,64,200,max=255), type="l", lwd=4)
for (t in 120:350){
	points(lowess(dataDim1$year,jitter(dataDim1$Dimension2), f=t/1000), col=rgb(255,64,64,15,max=255), type="l", lwd=.4)	
}
#
legend(1920,3.75,legend=c("1st Dimension","2nd Dimension"), col=c(rgb(30,144,255,200,max=255),rgb(255,64,64,200,max=255)), lwd=c(3,3), lty=c(4,1), bty="n")
#
axis(1, at = c(1920,1940,1960,1980,2000,2010), labels=c("1920","1940","1960","1980","2000","2010"))
axis(2, at = c(0:5), labels=c("0","1","2","3","4","5"), las=1)

############################################################################################################################################
# grey
## Plot
par(family="CMU Serif", mai=c(1.0,1,0.5,1.0))#,bg="grey95")
plot(dataDim1$year,dataDim1$Dimension1, ylim=c(0,4.9), xlim=c(1919,2010), type="l", lwd=2, lty=2, yaxt="none", xlab="Year", bty="n", col=rgb(255,64,64,1,max=255), xaxt="none", ylab="Number of Annual Initiatives")
rect(1919,-0.10,2012,6.4, col=rgb(247,247,247,150,max=255), border=NA)
points(jitter(dataDim1$year),dataDim1$Dimension2-0.04, col=rgb(140,140,140,40,max=255), pch=19)
points(jitter(dataDim1$year),dataDim1$Dimension1+0.04, col=rgb(140,140,140,40,max=255), pch=18, cex=1.4)
#
points(lowess(dataDim1$year,jitter(dataDim1$Dimension1), f=0.2), col=rgb(140,140,140,220,max=255), type="l", lwd=4, lty=3)
for (t in 120:350){
	points(lowess(dataDim1$year,jitter(dataDim1$Dimension1), f=t/1000), col=rgb(140,140,140,20,max=255), type="l", lwd=.4, lty=1)	
}
points(lowess(dataDim1$year,jitter(dataDim1$Dimension2), f=0.2), col=rgb(140,140,140,220,max=255), type="l", lwd=4)
for (t in 120:350){
	points(lowess(dataDim1$year,jitter(dataDim1$Dimension2), f=t/1000), col=rgb(140,140,140,20,max=255), type="l", lwd=.4)	
}
#
legend(1920,4.85,legend=c("1st Dimension","2nd Dimension"), col=c(rgb(140,140,140,220,max=255),rgb(140,140,140,220,max=255)), lwd=c(3,3), lty=c(3,1), bty="n")
#
axis(1, at = c(1920,1940,1960,1980,2000,2010), labels=c("1920","1940","1960","1980","2000","2010"))
axis(2, at = c(0:5), labels=c("0","1","2","3","4","5"), las=1)
points(1938,4.68,col=rgb(140,140,140,220,max=255),pch=18, cex=1.4)
points(1938,4.499,col=rgb(140,140,140,220,max=255),pch=19)









############################################################################################################################################
#			F	I	G	U	R	E 		3
############################################################################################################################################


############################################################################################################################################
# color
dataDiMM    <-  read.dta("Data_SPSR_2015_Leemann_orig.dta")
#
par(family="CMU Serif", mai=c(0.8,0.8,0.5,0.5))#,bg="grey95")
plot(dataDiMM$year,runmean(dataDiMM$political_elite ,k=35,alg="R",endrule="mean"), type="l", lwd=4, xlim=c(1920,2012),col=rgb(255,64,64,200,max=255), ylim=c(0,2.6),yaxt="none", xlab="Year", bty="n",  xaxt="none", ylab="Initiatives and Their Origins")
rect(1920,-0.10,2012,2.65, col=rgb(247,247,247,150,max=255), border=NA)
#
points(dataDiMM$year,dataDiMM$political_elite/1.9, col=rgb(255,64,64,20,max=255), pch=19)
points(dataDiMM$year,dataDiMM$civil_society/1.9, col=rgb(30,144,255,20,max=255), pch=19)
#
points(dataDiMM$year,runmean(dataDiMM$civil_society ,k=35,alg="R",endrule="mean"), type="l", lwd=4,col=rgb(30,144,255,200,max=255), lty=4)
points(dataDiMM$year,runmean(dataDiMM$political_elite ,k=35,alg="R",endrule="mean"), type="l", lwd=4,col=rgb(255,64,64,200,max=255))
	for (i in 1:400){
		points(dataDiMM$year,runmean(jitter(dataDiMM$political_elite , factor=2),k=35+(i-200)/40,alg="R",endrule="mean"), type="l",lwd=1, col=rgb(255,64,64,10,max=255))
		}
	for (i in 1:400){
		points(dataDiMM$year,runmean(jitter(dataDiMM$civil_society , factor=2),k=35+(i-200)/40,alg="R",endrule="mean"), type="l",lwd=1, col=rgb(30,144,255,10,max=255))
		}
#
axis(1, at = c(1920,1940,1960,1980,2000,2010), labels=c("1920","1940","1960","1980","2000","2010"))
axis(2, at = c(0,1/1.9,2/1.9,3/1.9,4/1.9,5/1.9), labels=c("0","1","2","3","4","5"), las=1)
#
legend(1920,2.5,legend=c("Parties and/or Politicians","Civil Society"), col=c(rgb(255,64,64,200,max=255),rgb(30,144,255,200,max=255)), lwd=c(5,5), lty=c(1,4), bty="n", cex=1.3)


############################################################################################################################################
# greyscale
par(family="CMU Serif", mai=c(0.8,0.8,0.5,0.5))#,bg="grey95")
plot(dataDiMM$year,runmean(dataDiMM$political_elite ,k=35,alg="R",endrule="mean"), type="l", lwd=4, xlim=c(1920,2012),col=rgb(140,140,140,200,max=255), ylim=c(0,2.6),yaxt="none", xlab="Year", bty="n",  xaxt="none", ylab="Initiatives and Their Origins")
rect(1920,-0.10,2012,2.75, col=rgb(247,247,247,150,max=255), border=NA)
#
points(dataDiMM$year,dataDiMM$political_elite/1.9+0.04, col=rgb(140,140,140,40,max=255), pch=19)
points(dataDiMM$year,dataDiMM$civil_society/1.9-0.04, col=rgb(140,140,140,40,max=255), pch=18,cex=1.4)
#
points(dataDiMM$year,runmean(dataDiMM$civil_society ,k=35,alg="R",endrule="mean"), type="l", lwd=4,col=rgb(140,140,140,200,max=255), lty=3)
points(dataDiMM$year,runmean(dataDiMM$political_elite ,k=35,alg="R",endrule="mean"), type="l", lwd=4,col=rgb(140,140,140,200,max=255))
	for (i in 1:400){
		points(dataDiMM$year,runmean(jitter(dataDiMM$political_elite , factor=2),k=35+(i-200)/40,alg="R",endrule="mean"), type="l",lwd=1, col=rgb(140,140,140,10,max=255))
		}
	for (i in 1:400){
		points(dataDiMM$year,runmean(jitter(dataDiMM$civil_society , factor=2),k=35+(i-200)/40,alg="R",endrule="mean"), type="l",lwd=1, col=rgb(140,140,140,10,max=255))
		}
#
axis(1, at = c(1920,1940,1960,1980,2000,2010), labels=c("1920","1940","1960","1980","2000","2010"))
axis(2, at = c(0,1/1.9,2/1.9,3/1.9,4/1.9,5/1.9), labels=c("0","1","2","3","4","5"), las=1)
#
legend(1920,2.5,legend=c("Parties and/or Politicians","Civil Society"), col=c(rgb(140,140,140,200,max=255),rgb(140,140,140,200,max=255)), lwd=c(5,5), lty=c(1,3), bty="n", cex=1.3)
points(1951,2.38,col=rgb(140,140,140,220,max=255),pch=19)
points(1939,2.27,col=rgb(140,140,140,220,max=255),pch=18, cex=1.4)









############################################################################################################################################
#			C	O	U	N	T	E	R	F	A	C	T	U	A	L	S
############################################################################################################################################



nbr10     <-  glm.nb(number_of_DD1 ~    change_pC + sharE + election + pre.election, data = mydata)

# Difference signatures: 1920 - 2011
x.prof.2011 <- data.frame(rbind(c(1.94412415,-1.0970883758,1,0),c(1.94412415,1.9781508522,1,0)))
names(x.prof.2011) <- c("change_pC", "sharE", "election", "pre.election")
predict(nbr10, x.prof.2011, type="response")

# Difference competition: 1920 - 2011
x.prof.2011 <- data.frame(rbind(c(min(change_pC),-1.0970883758,1,0),c(max(change_pC),-1.0970883758,1,0)))
names(x.prof.2011) <- c("change_pC", "sharE", "election", "pre.election")
predict(nbr10, x.prof.2011, type="response")






################################################################################################################ PLOT APPENDIX

par(family="CMU Serif", mfrow=c(2,1))
plot(c(1920:2011),attr(out.empty1, "prob.state")[,1],type="o", ylab=expression(paste("Pr(", S[t], "= k |", Y[t], ")")), xlab="Year", main="Model 1 Break: 2 Regimes", lwd=2, pch=19)
points(c(1920:2011),attr(out.empty1, "prob.state")[,2], col="blue",type="o", lwd=2, pch=19)
#points(c(1920:2011),attr(out.empty1, "prob.state")[,3], col="red",type="o", lwd=2, pch=19)
plot(c(1920:2011),attr(out.empty2, "prob.state")[,1],type="o", ylab=expression(paste("Pr(", S[t], "= k |", Y[t], ")")), xlab="Year", main="Model 2 Breaks: 3 Regimes", lwd=2, pch=19)
points(c(1920:2011),attr(out.empty2, "prob.state")[,2], col="blue",type="o", lwd=2, pch=19)
points(c(1920:2011),attr(out.empty2, "prob.state")[,3], col="green",type="o", lwd=2, pch=19)
#points(c(1920:2011),attr(out.empty2, "prob.state")[,4], col="red",type="o", lwd=2, pch=19)








########	
## lag explanatory variable

options(signif.symbols=c("***"=.01,"**"=.05,"*"=.1))
nbr1L <-  glm.nb(number_of_DD1[-c(1:4)] ~  change_pC[-c(89:92)] , data = mydata)
nbr2L <-  glm.nb(number_of_DD1[-c(1:4)] ~  sharE[-c(1:4)]  , data = mydata)
nbr3L <-  glm.nb(number_of_DD1[-c(1:4)] ~  change_pC[-c(89:92)] + sharE[-c(1:4)]  , data = mydata)
nbr4L     <-  glm.nb(number_of_DD1[-c(1:4)] ~   change_pC[-c(89:92)] + sharE[-c(1:4)] + election[-c(1:4)] , data = mydata)
nbr5L     <-  glm.nb(number_of_DD1[-c(1:4)] ~   change_pC[-c(89:92)] + sharE[-c(1:4)] + election[-c(1:4)] + pre.election[-c(1:4)] , data = mydata)


outi3 <- mtable("Model 1"=nbr1L,"Model 2" = nbr2L, "Model 3" = nbr3L, "Model 4" = nbr4L, "Model 5" = nbr5L, digits=2, summary.stats=c("McFadden R-sq.","N","Likelihood-ratio","p","BIC"))
outi3 <- relabel(outi3, "(Intercept)"="Constant", "change_pC"="Electoral Competition", "sharE"="Signature Requirements", "election"="Election Year", "pre.election"="Pre-Election Year")
toLatex(outi3)
outi3

# 10 year lag
options(signif.symbols=c("***"=.01,"**"=.05,"*"=.1))
nbr1L <-  glm.nb(number_of_DD1[-c(1:10)] ~  change_pC[-c(83:92)] , data = mydata)
nbr2L <-  glm.nb(number_of_DD1[-c(1:10)] ~  sharE[-c(1:10)]  , data = mydata)
nbr3L <-  glm.nb(number_of_DD1[-c(1:10)] ~  change_pC[-c(83:92)] + sharE[-c(1:10)]  , data = mydata)
nbr4L     <-  glm.nb(number_of_DD1[-c(1:10)] ~   change_pC[-c(83:92)] + sharE[-c(1:10)] + election[-c(1:10)] , data = mydata)
nbr5L     <-  glm.nb(number_of_DD1[-c(1:10)] ~   change_pC[-c(83:92)] + sharE[-c(1:10)] + election[-c(1:10)] + pre.election[-c(1:10)] , data = mydata)


outi3 <- mtable("Model 1"=nbr1L,"Model 2" = nbr2L, "Model 3" = nbr3L, "Model 4" = nbr4L, "Model 5" = nbr5L, digits=2, summary.stats=c("McFadden R-sq.","N","Likelihood-ratio","p","BIC"))
outi3 <- relabel(outi3, "(Intercept)"="Constant", "change_pC"="Electoral Competition", "sharE"="Signature Requirements", "election"="Election Year", "pre.election"="Pre-Election Year")
toLatex(outi3)
outi3




### Table with all votes

data.tab <- read.csv("Data_SPSR_2015_Leemann_Table.csv")
xtable(data.tab, digits=0)




##############################
# break model
##############################



# no breaks 
burnin=300000
mcmc <- 300000
out.emptyNEW0 <-  MCMCpoisson(number_of_DD1 ~ change_pC + sharE + election, data=mydata, m=0, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin, b0=0, B0=0.01, marginal.likelihood = "Laplace")  


# one break
burnin <- 300000
mcmc <- 300000
thin <- 300
startingval <- c(4.084e-01, 2.314e-01, 6.952e-02, -6.236e-05, 8.072e-01, 2.237e-01, -4.455e-01, 2.873e-01)
out.emptyNEW1 <-  MCMCpoissonChange(number_of_DD1 ~ change_pC + sharE + election, data=mydata, m=1, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin, marginal.likelihood = "none")#, beta.start=startingval)  

geweke.diag(out.emptyNEW1)
heidel.diag(out.emptyNEW1)
summary(out.emptyNEW1)
summary(out.emptyNEW1, quantiles=c(0.025,0.05,0.5,0.95,0.975))
plot(out.emptyNEW1)
attr(out.emptyNEW1, "prob.state")
plotState(out.emptyNEW1)
# even sharE_2 is significant on 90% ->>> sort(out.emptyNEW[,7])[3000*.95]

thin <- 1000
# two breaks 
burnin=300000
mcmc <- 300000
out.emptyNEW2 <-  MCMCpoissonChange(number_of_DD1 ~ change_pC + sharE + election, data=mydata, m=2, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin, marginal.likelihood = "Chib95")  

# three breaks 
burnin=100000
mcmc <- 200000
out.emptyNEW3 <-  MCMCpoissonChange(number_of_DD1 ~ change_pC + sharE + election, data=mydata, m=3, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin, marginal.likelihood = "Chib95")  

# four breaks 
burnin=100000
mcmc <- 200000
out.emptyNEW4 <-  MCMCpoissonChange(number_of_DD1 ~ change_pC + sharE + election, data=mydata, m=4, mcmc=mcmc, burnin=burnin, verbose=verbose,
                           thin=thin, marginal.likelihood = "Chib95")  



# Which is the best model
print(BayesFactor(out.emptyNEW0,out.emptyNEW1,out.emptyNEW2,out.emptyNEW3,out.emptyNEW4))



log.bayes <- BayesFactor(out.emptyNEW0,out.emptyNEW1,out.emptyNEW2,out.emptyNEW3,out.emptyNEW4)$BF.log.mat

colnames(log.bayes) <- c("No Break", "One Break", "Two Breaks", "Three Breaks", "Four Breaks")
rownames(log.bayes) <- c("No Break", "One Break", "Two Breaks", "Three Breaks", "Four Breaks")

xtable(log.bayes)



##############################
# plot
##############################


###### Cantonal data

dat1 <- read.dta("Data_SPSR_2015_Leemann_Kantone.dta")

mod1 <- lm( num_init ~ sig_init_rel, data=dat1[-1,])
mod2 <- glm( num_init ~ sig_init_rel, data=dat1[-1,],family="poisson")
mod3 <- glm.nb( num_init ~ sig_init_rel, data=dat1[-1,])


out.can <- mtable("Linear Regression"=mod1, "Poisson Regression" = mod2, "Negative binomial Model" = mod3, digits=2, summary.stats=c("R-squared","N","p","Log-likelihood","BIC"))
toLatex(out.can)


